#! /usr/bin/env python
#coding=utf-8
import math

class LinEqu: 
    def __init__(self, mat_aug):
        self.mat_aug = mat_aug;
	self.index = len(mat_aug[0])-1;
	self.solution = [0 for a in range(self.index)];
	print 'Original matrix: ';
	self.printMatrix();
    def eliminate(self):
	m = self.index;
	n = self.index+1;
	i,j=0,0;
	while True:
	    if i>m-2 or j>n-2: break;
	    print '\ni=%d, j=%d' % (i , j);
            mainElement = 0.0;
            for k in range(i+1, m):
                if math.fabs(mat_aug[k][i]>mainElement):
                    mainElement=mat_aug[k][i];
                    exchangeRow = k;
	    if not mainElement>0.0:
		print 'There is no solution or the solutions are more than 1. Exit.';
		return None;
                    
            for l in range(n):
		#Swap row i and exchangeRow that contains mainElement
                mat_aug[i][l], mat_aug[exchangeRow][l] = mat_aug[exchangeRow][l], mat_aug[i][l]; 
		mat_aug[i][l] /= mainElement;
            
            for u in range(i+1, m):
		pivot = mat_aug[u][j];
		for v in range(n): 
		    mat_aug[u][v]=mat_aug[u][v]-mat_aug[i][v]*pivot;

	    print 'After round %d of row swap and normalization: ' % i;
	    self.printMatrix();
	    i+=1;
	    j+=1;

    def backSub(self): 
	mat_aug = self.mat_aug;
	x=self.solution;
	m=self.index;
	n=self.index+1;
	b=[0 for row in range(self.index)];

	x[m-1]=mat_aug[m-1][n-1]/mat_aug[m-1][n-2];
	for i in reversed(range(self.index)):
	    b[i] = self.mat_aug[i][n-1];
	    sum = 0.0;
	    for j in reversed(range(i+1, n-1)):
		sum+=mat_aug[i][j]*x[j];
	    x[i] = (b[i]-sum)/mat_aug[i][i];

    def printMatrix(self):
	for row in range(len(self.mat_aug)):
	    print mat_aug[row];
	print;
    def printSolu(self):
	print
	print 'The soluiton is: ', self.solution;
    def main(self):
	self.eliminate();
	self.backSub();
	self.printSolu();

if __name__=='__main__':
    #mat_aug = [[9.0,3,4,7],[4,3,4,8],[1,1,1,3]];
    mat_aug = [[0.0,2.0,1.0,4.0],[1.0,1.0,2.0,6.0],[2.0,1.0,1.0,7]];
    equation = LinEqu(mat_aug);
    equation.main();
                
